-
Notifications
You must be signed in to change notification settings - Fork 76
Full implementation of AGN feedback for the cool-core destruction project #339
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
base: main
Are you sure you want to change the base?
Full implementation of AGN feedback for the cool-core destruction project #339
Conversation
|
@ShuangShuang0411 Some quick and general comments. I'll check the entire PR thoroughly after it is approved by other reviewers.
|
|
@jzuhone Since this PR is most related to your work, could you help review it? I'll also take a close look after you approve it. First, do you prefer to modify the existing |
|
@ShuangShuang0411 @hyschive I plan on reviewing this in detail over the next few days. |
example/test_problem/Hydro/ClusterMerger/yt_script/plot_slice-z.py
Outdated
Show resolved
Hide resolved
src/TestProblem/Hydro/ClusterMerger/Init_TestProb_ClusterMerger.cpp
Outdated
Show resolved
Hide resolved
src/TestProblem/Hydro/ClusterMerger/Flu_ResetByUser_ClusterMerger.cpp
Outdated
Show resolved
Hide resolved
src/TestProblem/Hydro/ClusterMerger/Init_TestProb_ClusterMerger.cpp
Outdated
Show resolved
Hide resolved
src/TestProblem/Hydro/ClusterMerger/Init_TestProb_ClusterMerger.cpp
Outdated
Show resolved
Hide resolved
src/TestProblem/Hydro/ClusterMerger/Par_Init_ByFunction_ClusterMerger.cpp
Outdated
Show resolved
Hide resolved
jzuhone
left a comment
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I left some comments.
src/TestProblem/Hydro/ClusterMerger/Par_Init_ByFunction_ClusterMerger.cpp
Outdated
Show resolved
Hide resolved
|
@ChunYen-Chen no, I think it's good for you to include those. |
|
@ChunYen-Chen I had a look at it and it looks good, so I merged. Anything else, aside from the conflict that I can fix? |
|
@jzuhone I think that's all. Thanks for the review. |
|
@hyschive I think this is ready for merge! |
|
@jzuhone @ChunYen-Chen Thanks for your great work updating and reviewing this PR! Please approve it if you think it’s ready. Since this is a fairly large PR, with 101 files modified, would you mind if I take another look before merging, at least for the files outside the |
| EoS_HTilde2Temp_CPUPtr, EoS_AuxArray_Flt, EoS_AuxArray_Int, | ||
| h_EoS_Table ); | ||
| # endif | ||
| if ( Temp <= 5e5 ) mass_cold[c] += fluid_acc[0]*dv; |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
@jzuhone Got it. If the value does not change, hard-coding should be fine.
hyschive
left a comment
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
This round of review address several compilation issues introduced by #458. I'll add more comments separately in the coming days.
src/TestProblem/Hydro/ClusterMerger/Flu_ResetByUser_ClusterMerger.cpp
Outdated
Show resolved
Hide resolved
src/TestProblem/Hydro/ClusterMerger/Flu_ResetByUser_ClusterMerger.cpp
Outdated
Show resolved
Hide resolved
src/TestProblem/Hydro/ClusterMerger/Flu_ResetByUser_ClusterMerger.cpp
Outdated
Show resolved
Hide resolved
src/TestProblem/Hydro/ClusterMerger/Flu_ResetByUser_ClusterMerger.cpp
Outdated
Show resolved
Hide resolved
src/TestProblem/Hydro/ExactCooling/Init_TestProb_ExactCooling.cpp
Outdated
Show resolved
Hide resolved
hyschive
left a comment
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
This round of review addresses additional compilation issues and includes a few minor comments.
hyschive
left a comment
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
This round of review focuses on the exact-cooling source files.
hyschive
left a comment
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
This round of review focuses on the ExactCooling test problem.
src/TestProblem/Hydro/ExactCooling/Init_TestProb_ExactCooling.cpp
Outdated
Show resolved
Hide resolved
src/TestProblem/Hydro/ExactCooling/Init_TestProb_ExactCooling.cpp
Outdated
Show resolved
Hide resolved
src/TestProblem/Hydro/ExactCooling/Init_TestProb_ExactCooling.cpp
Outdated
Show resolved
Hide resolved
…gamer-fork into cool-core-public * 'cool-core-public' of ssh://github.com/ShuangShuang0411/gamer-fork: Minor Update for CCD simulation Add more scripts Update for CCD simulations fix: prevent use of uninitialized variable chore(docs): Sync wiki to doc/wiki [skip-cd] Update output file
Co-authored-by: Hsi-Yu Schive <hyschive@gmail.com>
Co-authored-by: Hsi-Yu Schive <hyschive@gmail.com>
…gamer-fork into cool-core-public * 'cool-core-public' of ssh://github.com/ShuangShuang0411/gamer-fork: Apply suggestions from code review Apply suggestions from code review
* main: check pressure only for the EoSs that need check pressure only for the EoSs that need add a check for pressure in interpolation Default FLU_INT_SCHEME based on SUPPORT_SPECTRAL_INT
Co-authored-by: Hsi-Yu Schive <hyschive@gmail.com>
…r.cpp Co-authored-by: Hsi-Yu Schive <hyschive@gmail.com>
|
@hyschive I believe I have resolved all of your suggestions and I have merged |
hyschive
left a comment
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
This round of review focuses on the updated ClusterMerger test. Note that the review is not yet complete, and I'll add further comments on this test in the next few days.
| if ( FirstTime ) | ||
| { | ||
| const double dh_max = amr->dh[MAX_LEVEL]; | ||
| if ( R_acc/dh_max <= Threshold[0] ) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
| if ( R_acc/dh_max <= Threshold[0] ) | |
| if ( R_acc/dh_max <= Threshold[0] && MPI_Rank == 0 ) |
| // flag cells within the target radius, and if the radius is not resolved with a specific number (Threshold[0]) of cells | ||
| for (int c=0; c<Merger_Coll_NumHalos; c++) | ||
| { | ||
| if ( DIST_SQR_3D( Pos, CM_ClusterCen[c] ) <= SQR(25*R_acc) && R_acc/dh <= Threshold[0] ) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Should we regard 25 as a runtime parameter?
| { | ||
| const double dh_max = amr->dh[MAX_LEVEL]; | ||
| if ( R_acc/dh_max <= Threshold[0] ) | ||
| Aux_Message( stderr, "WARNING : MAX_LEVEL is less than the desired refinement level set in Input__Flag_User!! dh_max = %13.7e\n", dh_max ); |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
| Aux_Message( stderr, "WARNING : MAX_LEVEL is less than the desired refinement level set in Input__Flag_User!! dh_max = %13.7e\n", dh_max ); | |
| Aux_Message( stderr, "WARNING : MAX_LEVEL (%d) is less than the desired refinement level set in Input__Flag_User (R_acc = %13.7e, dh_max = %13.7e, Threshold[0] = %13.7e) !!\n", MAX_LEVEL, R_acc, dh_max, Threshold[0] ); |
| static char (*Merger_File_Prof)[1000]; // profile table of clusters | ||
| char (*Merger_File_Par) [1000]; // particle file of clusters |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I prefer replacing 1000 with MAX_STRING, which is defined in Macro.h and set to 512 by default.
| // --------------------------------------------------------------------------------------- | ||
| // (4) other variables | ||
| int AdjustCount = 0; // count the number of adjustments | ||
| int Merger_Coll_NumBHs; // number of BHs in the simulation |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
What's the difference between this variable and Merger_Coll_NumHalos? It's a bit confusing.
It also looks like Merger_Coll_NumHalos is mainly used during initialization (except for Flag_ClusterMerger.cpp). If so, should we replace it with Merger_Coll_NumBHs (or vice versa) for clarity?
| double (*CM_ClusterCen)[3]; // the center of each cluster | ||
| double (*CM_BH_Pos)[3]; // BH position of each cluster |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
What's the difference between CM_ClusterCen, CM_BH_Pos, and Merger_Coll_Pos? It's a bit confusing.
In the current implementation, it seems that CM_ClusterCen and CM_BH_Pos are always identical. However, I can imagine cases where the cluster center and the BH position might differ if their definitions change in the future. Adding some comments here would be helpful.
It also looks like Merger_Coll_Pos is only used during initialization. If so, should we replace it with CM_ClusterCen for clarity? Or at least we should add some comments (and perhaps rename Merger_Coll_Pos to for example Merger_Init_Pos for clarification)
| bool Flag = false; | ||
|
|
||
| // flag cells within the target radius, and if the radius is not resolved with a specific number (Threshold[0]) of cells | ||
| for (int c=0; c<Merger_Coll_NumHalos; c++) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Should it be Merger_Coll_NumHalos or Merger_Coll_NumBHs? I'm concerned because after a merger we only reduce Merger_Coll_NumBHs but not Merger_Coll_NumHalos .
| // Function : GetClusterCenter | ||
| // Description : Get the cluster centers | ||
| // | ||
| // Note : 1. Must enable Merger_Coll_LabelCenter |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Since we always call GetClusterCenter(), should we remove the Merger_Coll_LabelCenter option for clarity? Note that the code currently does not raise any error or warning when disabling Merger_Coll_LabelCenter.
| Aux_Message( stderr, "WARNING : file \"%s\" already exists !!\n", FileName ); | ||
|
|
||
| FILE *File_User = fopen( FileName, "a" ); | ||
| fprintf( File_User, "# ClusterCen_x/y/z : BH position\n" ); |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Clarify the units of all parameters.
hyschive
left a comment
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
This final round of review adds further comments on the updated ClusterMerger test. Once all comments are addressed, I will perform a final check before merge.
| // rho_gas : Average gas density within the accretion radius | ||
| // cs : Average sound speed within the accretion radius | ||
| // v_rel : Relative velocity of the gas within the accretion radius and the black hole |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
IIUC, these three variables refer to “hot” gas. If so, please clarify that.
| for (int j=0; j<PS1; j++) | ||
| for (int i=0; i<PS1; i++) | ||
| { | ||
| const double pos2[3] = { amr->patch[0][lv][PID]->EdgeL[0] + (0.5+i)*dh, |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Is there a particular reason for naming this variable pos2 rather than simply pos?
| if ( DIST_SQR_3D( pos2, CM_ClusterCen[c] ) > SQR(R_acc) ) continue; | ||
|
|
||
| double dr[3] = { pos2[0]-CM_ClusterCen[c][0], pos2[1]-CM_ClusterCen[c][1], pos2[2]-CM_ClusterCen[c][2] }; |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
DIST_SQR_3D already computes pos2[:]-CM_ClusterCen[c][:]. I suggest removing the redundant calculations.
| ang_mom[c][0] += dv * ( dr[1]*amr->patch[FluSg][lv][PID]->fluid[3][k][j][i] - dr[2]*amr->patch[FluSg][lv][PID]->fluid[2][k][j][i] ); | ||
| ang_mom[c][1] += dv * ( dr[2]*amr->patch[FluSg][lv][PID]->fluid[1][k][j][i] - dr[0]*amr->patch[FluSg][lv][PID]->fluid[3][k][j][i] ); | ||
| ang_mom[c][2] += dv * ( dr[0]*amr->patch[FluSg][lv][PID]->fluid[2][k][j][i] - dr[1]*amr->patch[FluSg][lv][PID]->fluid[1][k][j][i] ); |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
| ang_mom[c][0] += dv * ( dr[1]*amr->patch[FluSg][lv][PID]->fluid[3][k][j][i] - dr[2]*amr->patch[FluSg][lv][PID]->fluid[2][k][j][i] ); | |
| ang_mom[c][1] += dv * ( dr[2]*amr->patch[FluSg][lv][PID]->fluid[1][k][j][i] - dr[0]*amr->patch[FluSg][lv][PID]->fluid[3][k][j][i] ); | |
| ang_mom[c][2] += dv * ( dr[0]*amr->patch[FluSg][lv][PID]->fluid[2][k][j][i] - dr[1]*amr->patch[FluSg][lv][PID]->fluid[1][k][j][i] ); | |
| ang_mom[c][0] += dv * ( dr[1]*amr->patch[FluSg][lv][PID]->fluid[MOMZ][k][j][i] - dr[2]*amr->patch[FluSg][lv][PID]->fluid[MOMY][k][j][i] ); | |
| ang_mom[c][1] += dv * ( dr[2]*amr->patch[FluSg][lv][PID]->fluid[MOMX][k][j][i] - dr[0]*amr->patch[FluSg][lv][PID]->fluid[MOMZ][k][j][i] ); | |
| ang_mom[c][2] += dv * ( dr[0]*amr->patch[FluSg][lv][PID]->fluid[MOMY][k][j][i] - dr[1]*amr->patch[FluSg][lv][PID]->fluid[MOMX][k][j][i] ); |
|
|
||
| for (int c=0; c<Merger_Coll_NumBHs; c++) | ||
| { | ||
| MPI_Allreduce( ang_mom[c], ang_mom_sum[c], 3, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD ); |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
It seems unnecessary to declare ang_mom_sum as a global variable in Init_TestProb_ClusterMerger.cpp. How about making it a local variable in this function instead?
| real P_old_sqr = SQR(fluid[MOMX]) + SQR(fluid[MOMY]) + SQR(fluid[MOMZ]); | ||
| real Ekin_old = ( dens_old == 0.0 ) ? (real)0.0 : 0.5*P_old_sqr/dens_old; | ||
| real P_new = SQRT( 2*fluid[DENS]*(EngySin + Ekin_old) ); | ||
| double JetSign = SIGN( Vec_c2m[0]*CM_Jet_Vec[c][0] + Vec_c2m[1]*CM_Jet_Vec[c][1] + Vec_c2m[2]*CM_Jet_Vec[c][2] ); |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
| real P_old_sqr = SQR(fluid[MOMX]) + SQR(fluid[MOMY]) + SQR(fluid[MOMZ]); | |
| real Ekin_old = ( dens_old == 0.0 ) ? (real)0.0 : 0.5*P_old_sqr/dens_old; | |
| real P_new = SQRT( 2*fluid[DENS]*(EngySin + Ekin_old) ); | |
| double JetSign = SIGN( Vec_c2m[0]*CM_Jet_Vec[c][0] + Vec_c2m[1]*CM_Jet_Vec[c][1] + Vec_c2m[2]*CM_Jet_Vec[c][2] ); | |
| const real P_old_sqr = SQR(fluid[MOMX]) + SQR(fluid[MOMY]) + SQR(fluid[MOMZ]); | |
| const real Ekin_old = ( dens_old == (real)0.0 ) ? (real)0.0 : (real)0.5*P_old_sqr/dens_old; | |
| const real P_new = SQRT( (real)2.0*fluid[DENS]*(EngySin + Ekin_old) ); | |
| const real JetSign = SIGN( Vec_c2m[0]*CM_Jet_Vec[c][0] + Vec_c2m[1]*CM_Jet_Vec[c][1] + Vec_c2m[2]*CM_Jet_Vec[c][2] ); |
| real P_new = SQRT( 2*fluid[DENS]*(EngySin + Ekin_old) ); | ||
| double JetSign = SIGN( Vec_c2m[0]*CM_Jet_Vec[c][0] + Vec_c2m[1]*CM_Jet_Vec[c][1] + Vec_c2m[2]*CM_Jet_Vec[c][2] ); | ||
|
|
||
| double P_old_perp[3], P_old_para, P_new_para; |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
| double P_old_perp[3], P_old_para, P_new_para; | |
| real P_old_perp[3], P_old_para, P_new_para; |
| fluid[MOMY] += CM_BH_Vel[c][1] * fluid[DENS]; | ||
| fluid[MOMZ] += CM_BH_Vel[c][2] * fluid[DENS]; | ||
|
|
||
| fluid[ENGY] = 0.5 * (SQR(fluid[MOMX]) + SQR(fluid[MOMY]) + SQR(fluid[MOMZ])) / fluid[DENS] + |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
| fluid[ENGY] = 0.5 * (SQR(fluid[MOMX]) + SQR(fluid[MOMY]) + SQR(fluid[MOMZ])) / fluid[DENS] + | |
| fluid[ENGY] = (real)0.5 * (SQR(fluid[MOMX]) + SQR(fluid[MOMY]) + SQR(fluid[MOMZ])) / fluid[DENS] + |
|
|
||
| Cool-core destruction: | ||
| ======================================== | ||
| Reproduce the cool-core destruction results. Shuang-Shuang Chen (arXiv:2412.13595) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
| Reproduce the cool-core destruction results. Shuang-Shuang Chen (arXiv:2412.13595) | |
| Reproduce the cool-core destruction results of Chen et al. 2026 (arXiv:2412.13595) |
| 0. Please contact Prof. Hsi-Yu Schive regarding the unreleased exact cooling feature. | ||
| 1. `sh download_ic_cool_core_destruction.sh` | ||
| 2. `sh setup_input_files.sh <case_number>` | ||
| --> For the high resolution case, please link or copy the input files from |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
| --> For the high resolution case, please link or copy the input files from | |
| --> For the high-resolution case, please link or copy the input files from |
This pull request includes the complete implementation of AGN feedback for the cool-core destruction project. The key components are:
ClusterMergertest problem.ExactCoolingtest problem for the currently privateExactCoolingsource term module.